% March, 2012
%
% Computes the scattering of a plane wave with two objects
% kh    Helmholtz parameter
% angle Angle of the incoming wave

function twoobj(kh,theta_in)

check_sol = true;
verif2d   = true;

% number of points per wavelength
ppl = 20;
acc = 1e-10;

nobj            = 2;
% The number of points on each contour.
n               = 200;  
ntot            = n*nobj;

% Grid size for 2d visualization
ngrid           = 250;

contours        = cell(nobj,1);
Ctot            = [];
hmax_phy        = 0;

% First contour
[C,len,xin,tmp] = get_geometry_star(n,[0;0],5,0.);
contours{1}     = C;
h_phy           = max(sqrt(C(2,:).^2 + C(5,:).^2));    
hmax_phy        = max(hmax_phy,h_phy);
Ctot            = [Ctot,C];          

% Second contour
[C,len,xin,tmp] = get_geometry_star(n,[4;0],5,0.);
contours{2}     = C;
h_phy           = max(sqrt(C(2,:).^2 + C(5,:).^2));    
hmax_phy        = max(hmax_phy,h_phy);
Ctot            = [Ctot,C];          

xctot       	= mean(Ctot(1,:));
yctot       	= mean(Ctot(4,:));
Rtot            = max(sqrt((Ctot(1,:) - xctot).^2 + (Ctot(4,:) - yctot).^2 ));

lambda = 2*pi/kh;

if verif2d
  figure(2), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
  figure(3), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
end

% Check the size of the discretization vs the wavelength
fprintf(1,'h / lambda = %g\n', kh*hmax_phy/(2*pi))

if kh*hmax_phy/(2*pi) > 1./10
  error('mesh too coarse')
end

if verif2d
  xx      = linspace(xctot - 3.5*Rtot,xctot + 3.5*Rtot,ngrid);
  yy      = linspace(yctot - 3.5*Rtot,yctot + 3.5*Rtot,ngrid);
  [xxx,yyy] = meshgrid(xx,yy);
  %%% Determine which nodes in the visualization mesh are inside the contour
  %%% by evaluating Cauchy integrals.
  %%% Note that one integral is evaluated for each node in the mesh so this
  %%% can be slightly expensive. 
  %%% If you want a finely resolved visualization mesh, it might be worth
  %%% it to precompute "ind" and store it.
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];
  nz  = size(zzz,2);
  zz =     (ones(nz,1)*Ctot(1,:) - zzz(1,:)' * ones(1,ntot)) + ...
	   1i * (ones(nz,1)*Ctot(4,:) - zzz(2,:)' * ones(1,ntot));
  h   = 2*pi / n;
  dz  = h*Ctot(2,:)' + h*1i*(Ctot(5,:)');
  uu  = (1/(2*pi*1i))*(1./zz)*dz;
  ind = find(abs(uu) > 0.01);
else
  % Verify the solution on a circle
  theta_zzz = linspace(0,2*pi,500).';
  xxx       = xctot + 3*Rtot*cos(theta_zzz);
  yyy       = yctot + 3*Rtot*sin(theta_zzz);
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];

  figure(1)
  plot(zzz(1,:),zzz(2,:),'k')		
end


% Assemble the matrix
matsA    = cell(nobj,1);
matA   = zeros(ntot);
for i=1:nobj
  i0 = (i-1)*n + 1; i1 = i0 + n - 1;
  matAii = get_A_single_diag(contours{i},1:n,kh);
  matsA{i}          = matAii;
  if (check_sol)
    for j=1:nobj
      i0 = (i-1)*n + 1; i1 = i0 + n - 1;
      j0 = (j-1)*n + 1; j1 = j0 + n - 1;
      if (i == j)
	matA(i0:i1,j0:j1) = matAii; 
      else
	matA(i0:i1,j0:j1) = get_A_offd_noquad(Ctot,i0:i1,j0:j1,kh);
      end
    end
  end
end

[contours_up,matsA_up,nobjsx_new,nobjsy_new,Dplt] = factor_level(2,1,contours,matsA,kh,acc,ppl,2);


% Direct solution
v     = -incoming_wave(kh,theta_in,Ctot(1,:),Ctot(4,:)); v = v.';
sigma = matA\v;
phi_scattered_ref  = evalpot(zzz,Ctot,sigma,kh);

% Approximate solution
v     = -incoming_wave(kh,theta_in,contours_up{1}(1,:),contours_up{1}(4,:)); v = v.';
sigma = matsA_up{1}\v;
phi_scattered      = evalpot(zzz,contours_up{1},sigma,kh);

err = phi_scattered - phi_scattered_ref;

if verif2d
  step     = 4;
%  [Inew,idx] = sort(Is);
 % Cplot      = Cs(:,idx);
%  Cplot      = Cplot(:,1:step:end);
     
  figure(2), hold all;
  err(ind) = NaN;
  err      = reshape(err,ngrid,ngrid);
  contourf(xxx,yyy,log10(abs(err)),-16:2:0)
  cblabel = colorbar('vert');
%  plot(Cplot(1,:),Cplot(4,:),'b.','markersize',20)
%  plot(D(1,:),D(2,:),'r','linew',2)
  axis equal
  axis([-10 10 -10 10])
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);


  plot(Dplt{1}(1,:),Dplt{1}(2,:),'r.');
  plot(Dplt{2}(1,:),Dplt{2}(2,:),'r.');

  % set(tlabel,'FontName', 'cmr10','FontSize',15)
%  print(gcf, '-dpng', '2obj_err.png');
  matlab2tikz('2objerr.tex')
  
  phi_in		= incoming_wave(kh,theta_in,xxx,yyy);
  phi_scattered_ref	= reshape(phi_scattered_ref,ngrid,ngrid);
  phi_tot		= phi_in + phi_scattered_ref;
  phi_tot(ind)		= NaN;
  figure(3), hold all;
  contourf(xxx,yyy,real(phi_tot))
  cblabel = colorbar('vert');
  % plot(Cplot(1,:),Cplot(4,:),'r.','markersize',20)
  axis equal
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);
  % set(tlabel,'FontName', 'cmr10','FontSize',15)
  %  print(gcf, '-dpng', '2obj_phitot.png');
%  matlab2tikz('2objphitot.tex')
else
  fprintf(1,'\n\n ERROR :  %g \n',max(abs(err)))
end

figure(4), hold all;

[D,h_D]          = get_circle(contours_up,1,kh,ppl,1.);
[Is,Cs,V,matP,D] = get_skeleton(kh,contours_up{1},matsA_up{1},acc,D,h_D);

% Approximate solution
v     = -incoming_wave(kh,theta_in,Cs(1,:),Cs(4,:)); v = v.';
sigma = matP*v;
phi_scattered      = evalpot(zzz,Cs,sigma,kh);
phi_scattered      = reshape(phi_scattered,ngrid,ngrid);

err = phi_scattered - phi_scattered_ref;

if verif2d
  step     = 4;
  [Inew,idx] = sort(Is);
  Cplot      = Cs(:,idx);
  Cplot      = Cplot(:,1:step:end);
  
  figure(4), hold all;
  err(ind) = NaN;
  contourf(xxx,yyy,log10(abs(err)),-16:2:0)
  cblabel = colorbar('vert');

 % plot(Cplot(1,:),Cplot(4,:),'b.','markersize',20)
  plot(D(1,:),D(2,:),'r.')
  axis equal
  axis([-10 10 -10 10])
  view(0,90)
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);
  % set(tlabel,'FontName', 'cmr10','FontSize',15)
 % print(gcf, '-dpng', '2obj_err2.png');
  matlab2tikz('2objerr2.tex');
else
  fprintf(1,'\n\n ERROR :  %g \n',max(abs(err)))
end
keyboard
return
